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Abstract 

Let i? be a discrete valuation ring (DVR) and K be its fraction field. If M is a matrix over R admitting 
a LU decomposition, it could happen that the entries of the factors L and U do not lie in R, but just in K. 
^— ■») , Having a good control on the valuations of these entries is very important for algorithmic applications. In 

■ the paper, we prove that in average these valuations are not too large and explain how one can apply this 
result to provide an efficient algorithm computing a basis of a coherent sheaf over from the knowledge 

^ , of its stalks. 
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Throughout the paper, we fix a ring R equipped with a discrete valuation vji : i? — > N U {oo}. We 
PsJ . assume that vr is normalized so that it takes the value 1 and we fix an element vr e i? such that vb,{t^) = 1. 

We also assume that R is complete with respect to the distance defined by vu- The residue field of R and 
its fraction field are denoted by k and K respectively. The valuation vr extends uniquely to K and we 
, , continue to denote this extension by vr. We finally set q = Card k and assume that q is finite. Two typical 

r> . examples of this are (1) i? = Zp (the ring of p-adic integers) equipped with the usual p-adic valuation and 

■ (2) i? = k[[x]] (the ring of power series) where A: is a finite field. 

If c? is a positive integer, we denote by il the ring of square matrices of size d with coefficients in R. 
It is a compact additive group whose Haar measure is denoted by ^. We assume that /i is normalized so 
that is a probability space (i.e. fi{il) = 1). Thus, it makes sense to study some statistics on O. 

Surprinsingly, the literature around this subject seems to be very poor. Nevertheless related questions were 
already addressed by Abdel-Ghaffar in [ 1 1 and Evans in |4|: the main result of 1 1] is the computation of the 
law of the random variable "valuation of the determinant" (in the case where i? is a power series ring but 
his argument works for a more general discrete valuation ring) whereas, in [4 |, Evans studies the random 
variable "valuation of the elementary divisors". 

In this paper, we are mainly interested in the random variable Vl- "valuation of the L-part in the 
LU decomposition". We give several estimations of its law, its expected value and its standard deviation. 
Roughly speaking, we prove that K[Vl] = logg d + 0{l) and (t{Vl) = 0(1) where the notation 0(1) refers 
to a quantity bounded by a universal constant. We also bound from above the probability that Vl deviates 
from its expectation. For more precise statements, we refer to Theorem ll.il Theorem ll.2l and Corollarv ll.3l 
in the introduction of SjT] 

In Sj2] we move to algorithmic applications. Firstable, we propose in i j2.1l a stable algorithm to compute 
a LU decomposition of a matrix over R (unfortunately standard Gauss elimination is far for being stable) 
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and analyze closely the losses of precision it generates in average (which turn out to be optimal in some 
sense). i j2.2l is devoted to a study of the notion of "simultaneous PLU decomposition", which will play an 
important role for our next (and main) application presented in m.3\ This application is of geometric nature. 
We let X = A]^ be the affine line over K. Recall that a coherent subsheaf of 7^ C Oj^ (where d is some 
integer) is determined by the data of all its stalks J'x C Oj^ ^ at all closed points x E X. Furthermore, we 
know that all such sheaves T as above admit a global basis. In m.3\ we describe an algorithm that computes 
a basis of J' knowing all its stalks and, once again, analyze its stability (which will turn out to be rather 
good). 

1 Some statistics related to LU decomposition 

If A is a (commutative) ring, we shall denote by Md (A) the ring of square dxd matrices. Recall that we have 
endowed fi = Md{R) with its Haar measure. Choosing a matrix at random with respect to this measure 
is just choosing independently each entry at random with respect to the Haar measure on R. Furthermore, 
since R is complete, every element x G R can be written uniquely as an infinite sum x = X^i^o '^^'^^ where 
the coefficients a^'s are taken in a fixed set 72. C i? of representatives of elements of k (i.e. the restriction to 
TZ of the canonical projection i? — > fc is bijective) and conversely, any sum X^i^o '^i'^^ above converges 
and then defines an element in R. With this description, generating a random element (with respect to the 
Haar measure) of R is just choosing at random all coefficients a/s in TZ independently and uniformly. 

We shall say that a matrix M G Md{K) admits a LU decomposition if it can be factorized as a product 
L{M) ■ U{M) where: 

• L{M) is a uni{3 lower triangular matrix with coefficients in K, and 

• U{M) is a upper triangular matrix with coefficients in K. 

We underline that, even if M has coefficients in R, we do not require that L{M) and U{M) belong to 
Md{R). Here are some well known facts: (1) an invertible matrix AI admits a LU decomposition if and 
only if all its principal minors do not vanish and (2) when it exists, a LU decomposition is unique {i.e. the 
matrices L{M) and U{M) are uniquely determined). We will consider L and U as two partially defined 
functions on Md{K). For lu E such that L{uj) is defined, let us denote by Vl(i^) the opposite of the 
smallest valuation of an entry of L{uj). The aim of this section is to study the random variable Vl- Here are 
the main results we will prove. 

Theorem 1.1. Setting 

oo 

= (1) 

v=l 

we have E{q, d) - < ¥\Vl\ ^ E{q, d). 

Furthermore, the distance between E{q, d) and log^ d is bounded by ^^^(^2) ('^^d by 1 if q ^ 3j. 

Theorem 1.2. For all positive real number I, we have: 

P[|Fi-l0g^d-i| >^+i] ^ (2 + Mogq). 

Corollary 1.3. The standard deviation o/Vl is bounded by an explicit universal constant (which can be 
chosen equal to 6.5). 

1.1 Some useful tools 

This subsection gathers some preliminaries to the proof of Theorem IL II Theorem IL2I and Corollarv ll.3l 
We first recall some basic facts about LU decomposition, then introduce the random variables Vij's (which 
will play a crucial role in the sequel) and finally prove several important properties of them. 

'it means that all diagonal entries are equal to 1. 
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1.1.1 Cramer's rule for LU decomposition 



Let M e Md{K). A useful formula for our purpose is an analogue of Cramer's rule which gives a closed 
expression of the entries of L{M) as a quotient of two determinants. This formula appears for instance in 
||6|, §1.4; let us recall it briefly. If / and J are two subsets of {1, . . . , d}, we denote by Af/.j the submatrix 
of M obtained by deleting all columns and rows whose index are not in / and J respectively. The i-th 
principal minor is then the determinant of the matrix Mjj where / = {I, . . . we will denote it by 
Si(M). With these notations, we have: 

ifi>j : LiM'U^j ^ '^T{M) where /={!,..., and J = {!,..., j} (2) 

iii^j : U{M\^, = f^^';' where/ = {l,...,i}andJ = {!,..., I- (3) 

The proof of these formulas is not difficult. For Formula (|2]i, note that L(M)/.jC/(M)jj — A//, j provided 
that J has the particular shape J = {!,... then, passing to the determinant, we get dct i(A/)/_j • 
det U{M)j^j = det Mj j and the desired relation follows by combining these equaUties for I = J and 
I = {1, . . . , j — l,i}. The proof of Formula (O is similar. 



1.1.2 The random variables Vi, j 

The aim of this paragraph is to define a collection of mututally independent random variables Vij : fl 
NU{cx)}(l ^ j ^ d); they will be very useful in the sequel to study Vl- We first construct a collection 
of random variables Xij : ^ R {1 ^ i ^ j ^ d). The construction goes by induction on j. We start with 
a matrix ui in fi. We first define Xi^i{uj) to be the top left entry of oj. We then enter in the second round (i.e. 
j — 2). As before, we begin by letting Xi^2{^^) denote the (1, 2)-th entry of cj but, before defining W2,2(w) 
we do the two following modifications on the matrix lj: 

• if the valuation of 2(w) is less than the valuation of Xi we swap the two first columns of oj 
and, tiien 

• we clear the (1, 2)-th entry of oj by adding to its second column a suitable multiple of its first column 
(note that it is always possible because if the top left entry — which serves as pivot — vanishes, so 
does the (1, 2)-th entry). 

Doing these operations, the coefficient of in position (2, 2) may have changed and we define X2,2(<^) to 
be the new (2, 2)-th entry of u. The general induction step works along the same ideas. Assume that, after 
the {j — l)-th step, we have ended up with a matrix uj such that un'j' = when i' < j' < j. We define 
Xi,j (w) by induction on i by applying the following process successively for i = 1, 2, . . . , j — 1: 

• first, we set Xij (oj) to the {i, j)-th entry of (the current) oj; 

• second, if the valuation of Xij{oj) is less than the valuation of the {i, i)-th entry of (the current) w, 
we swap the first row of oj with its i-th one; 

• third, we clear the {i,j)-th entry of oj by adding to its j-th column a suitable multiple of its first 
column. 

We finally let Xj,j{uj) denote the j-th diagonal entry of (the current) oj. For all with 1 ^ i ^ j ^ d, 
we also set Vij = vji{Xi,j) and Vi = Vi^i. The V^./s take values in N U {oo} and they are finite almost 
everywhere. Algorithm[T] summarizes the construction of the V^j 's. 

Proposition 1.4. The random variables Xij (I ^ i ^ j ^ d) are uniformely distributed and mutually 
independent. 

Proof. Set I — { {i , j) \1 ^ i ^ j ^ d} . Suppose we are given a family x = of elements of 

R. We consider the following set: 

n{x) = {oj en\ Xij{oj) = x^^j, e / }. 
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Algorithm 1: The construction of the random variables Vij's 

Notation: : * Wij denotes the (i, j )-th entry of a; 
: ★ ujj denotes the j-th row of uj 

1 for j from Ito ddo 

2 for ifrom 1 fo j — 1 do 

4 T£vii{uji^j) < vji{uji,i) then swap uij and uji; 

5 it 7^ then ujj ojj — • uji; 

6 Vj^j ^ VRiujjj); 



Set — vji{xij) and for all i, let ji denote the first index in {1, . . . , d + 1 — i} such that Vij-i+j^ is 
equal to mm{vi,i, . • . , w^.d)- This sequence of integers (ji) is the code of a certain permutation a of 

{1, . . . ,d} defined by the following rule. We write all the integers between 1 and d. We define cr(l) to be the 
ji-th written integer (that is ji) and we erase it. We then define a{2) to be the j2-th integer which remains 
written (that is j2 is j2 < ji and j2 + 1 otherwise), we erase it and we continue. Let denote the subset of 
{1, . . . , d}^ consisting of couples such that i > a~^{j). One can check that it has cardinality 
Consider the function : il{x) — > mapping uj to the family ■ Following the construction 

of the j's, one can check that fx is a bijection. 

Now, we globalize the previous construction. Let C/ be a subset of containing a distinguished element 
X and such that ly = for all y ^ U. With this assumption, the collection of functions fy's {y varying in 
U) defines a bijection between fl{U) = {oj E n \ ^j^/ & U} and U x i?^^. It is morever easy 

to check that this bijection preserves the measure; in other words 

P[(X,,)(,,,)e/eC/]=M(f/) (4) 

where /i denotes the Haar measure on Rf But, since the function vr is locally constant on i?\{0}, any 
open subset U C (i?\{0})^ can be written as a disjoint union of subsets U' on which the function y ^ ly 
is constant. Therefore the equality ^ holds for all these U . Since furthemore the complement of (i?\{0})^ 
in R^ is a measure-zero set, the equality (HI holds for all open subset of i?^. □ 

Corollary 1.5. The random variables Vij (1 ^ i ^ j ^ d) are mutually independent and they all follow a 
geometric law of parameter 1 — q^^ (i.e. they take value v with probability (1 — q)q^~^). 

Proof. Clear after Proposition IL4I □ 

Another interest of the V^ j's is that they are closely related to Vl- The following Proposition precises 
this relationship. 

Proposition 1.6. We have ma.x(Vi, V2, ■ ■ ■ , Vd) — VR{det) ^ Vl ^ max(Vi, . . . , Vd-i) (recall that Vi = 
Vi^i by definition). 

Proof. Let w G il. To avoid confusion, agree to call Tj (w) the matrix w computed by Algorithm[T|(run with 
UJ as input) after the j-th iteration of the main loop and reserv the notation lo for the matrix we have started 
with. It follows from the construction that Tj{uj) has the following particular shape: if i' < j' ^ j, then 
the (z', j')-th entry of cuj vanishes. Moreover, clearly, Tj{uj) is obtained from uj by performing successive 
elementary operations on the first j columns. Therefore, if J = {1, . . . , j} and if / is a subset of {1, . . . , d} 
of cardinality J, we have detcj/ j = ± det(wj )/,,/. In particular these two determinants have the same 
valuation. Fix a couple [i,]) such that 1 ^ j ^ z ^ d and set / = {1, . . . , j — 1, z}, J = {1, . . . , j}. From 
Formula (|2]i and what we have said before, we derive: 

VR{L{oj)tj) = vii{detTj{u!)ij)-VR{detTj{u!)j,j) 

= VR{Tj{uj)ij) - VR{Tj{uj)j,j) = VR{Tj{uj)ij) - Vj{uj). 

Since all coefficients of ujj lie in R, so does its determinant. It follows that VR{Tj{u!).ij) ^ and conse- 
quently that VR{L{uj)ij) ^ —Vj{uj), which proves the second inequality. To establish the first one, note 
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that Lo and Tj{uj) share the same determinant up to a sign. Thus there must exist an index i, necessar- 
ily not less than j, such that V]^{Tj{uj)i,j) ^ VR{detLu). For this particular i, we have vii{L{uj)ij) ^ 
i'it:(det oj) — and then Vl(cli) ^ t/, (oj) — Ufl(det oj). The conclusion follows. □ 

Remark 1.7. In the same way, we can prove that the valuation of the i-th minor of a; G is equal to 
X]i=i ■ • • 1 Combining this with Corollarv 11.51 one can easily recover 

Abdel-Ghaffar's formula X^iLi ^t3T (see Theorem 3 of UJ) giving the expected value of the random vari- 
able "valuation of the determinant". 

1.2 Proof of the main results 
1.2.1 Estimation of the expected value 

This subsection is devoted to the proof of Theorem ll.il 



Estimation of the expected value of Vl Let V = max(Vi, V2, . . . , Vd). The event "V < w" occurs if 
and only if Vi,i < v for all index i, and Corollarv II .51 shows that it happens with probability (1 — q~'"Y. 
The expected value of V is then equal to X^^i v] — [l ^ (1 ^ 9 that is exactly E{q^ v). 

On the other hand. Proposition [131 implies that E[F] - E[u/? (det)] ^ ¥.\Vl\ ^ '&\V]. Moreover, by Abdel- 
Ghaffar's Theorem, we know that the expected value of w;j(det) is given by 5*31 hence is less 
that < YlTLi = The first part of Theorem ll.ll is proved. 



q- 



Estimation of E{q, v) Consider the function / :xi~->l — (1 — g ^) . It is decreasing on the interval 
[0, 00) and therefore one can write: 



00 POO 



f{x)dx ^ E{q,d) ^ / f{x)dx ^ -1+ f{x)dx. 
Ji Jo 



Doing the substitution y = 1 — q ^, we get: 







f{^)dx = r— • / dy^-—- {l+y + y' + ... + y'^-')dy 



logg Jo l-y logq J a log 9 



where Hd ~ l + i + + the harmonic series. It is well known that 7 + logd ^ Hd ^ 1 + logd 
where 7 is the Euler's constant. Therefore E{q, d) is almost equal to log^ d, the error term being bounded 
by a universal constant. The second part of Theorem ll.ll follows. 



Some additional remarks We would like first to emphasize that the difference E{q, d) — log^ d does not 
converge to when q and/or d goes to infinity. Indeed the following Lemma shows that, when log^ d is far 
from an integer and q is large, E{q, d) might be closer to the integral part of log^ d than to log^ d itself. 

Lemma 1.8. For all q qnd d, 

\E{q,d)-[\0g^d\\<-^-q-'^<'°^^''^'''> 

where [log^ d] and dist(log^ d, N) denotes respectively the integral part and the distance to N o/log^ d. 
Proof. We claim that the function / : a; 1 — (1 — q^'-^ Y satisfies: 

1 - ^ < f{x) < for all x^Q. (5) 

d ^ ' q"" 

Indeed, the second inequality directly comes from the standard inequality + ^ 1+td whereas the first 
one is a consequence of AM-GM inequality applied with the numbers dq^^ and 1 — g^^, 1 — g"^, . . . , 1 — 
q-"" (d times). If vq = [log^ d], we then get vq + Y.7=i ^ < E{q, d) ^ vq + Y.T=vo+i ^' '^^i'^^ g^^^^- 

q q"^ r., q d 



E{q,d)-vo s$ 



q-1 d ^ ' q-1 q^o+i 

The Lemma follows from this. □ 
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Let us end this paragraph by a last remark: the sum E{q, d) can also be exactly computed. Indeed, we 
have: 

■u=l v = lk=l ^ ^ fe=l ^ ^ 

Nevertheless, this expression does not yield the order of magnitude of E{q,d); indeed, each term in the 
latter sum (the one over k) can individualy be very large whereas the sum itself grows rather slowly. 



I. 2.2 Estimation of the law of Vl 

We now start the proof of Theorem ll.2l The strategy is quite clear: we use Corollary [T3] and Proposition 

I I. 61 to bound from below and from above the distribution function of Vl- First, let us investigate the con- 
sequences of the inequaUty Vl ^ V (where we recall that we have set V = max(Vi, . . . , Vd))- For all 
(nonnegative) real number x, it implies that: 

d 

P[Vl <x]^ P[V < x] = Y[V[V <x]^{l- q-^'f ^1-d- q-\ (6) 

1=1 

It is a bit more tricky to use the other inequality Vl ^ V — VLt{dei) because Ufl(det) and the V^'s are 
certainly not independent (c/ Remark fLTl i. Nevertheless, one can pick two nonnegative real numbers x and 
t and consider the event E^^t ■ > x + t and (det) ^ t". It is clear that V — Vfi {det) is always greater 
than X when E^^t occurs. Thus we have: 

P[Vl ^x]i^ V[V - Ufl(dct) ^ x] P[E^,t] i^l-P[V i^x + t]- ¥[vR{det) > t]. (7) 

Moreover we know that P[V ^ x + t] ^ (1 — q-x-tyi from Abdel-Ghaffar' s result (see fl]), we 

-t + 2 

derive F[Vd > t] ^ "^q-i ■ Indeed, Abdel-Ghaffar Theorem states that for all integer v, the equality 
P[-i;fl(det) ^ -y] = (1 - q^"^^){l - q"""^) • • • (1 - q-'"-'^) holds. In particular P[wfl(det) < ^ 
1 — Yli=i 9"" * ^ 1 — ^Tij- Taking v — [t], we get the claimed result. Putting these inputs in (|7]i, we 

P[Vl ==; a;] s; 1 - (1 - q-"^-*)'' ^ 



^d 

obtain: 



' 9-1 

This estimation being true for all t, one can optimize it on t. For simplicity, let us define u = 1 — g"^"*; 
the variable u now varies in [1 — q^^, 1], and for all u in this range, one have P[Vl > v] ^ 1 — f{u) where 
f{u) ~ u + dX{l — u), X — /(q-i) ■ Assume that A < 1. A quick study of / shows that it is minimal when 
u ^ uq = A^/ (''^i) . Moreover, one can check (using AG-MG inequality for instance) that uq always lies in 
the interval [1 - g"^, 1]. It follows that P[Vl < x] ^ 1 - /(wq) = 1 - A • (d - (d - l)uo). We can further 
simplify this formula and write a bound depending only on A. For this, remark that A ^ (1 + ^fr^)''^- 
Raising to the power d — 1, we find uq ^ 1 + ^n*! ^^n: 

F[Vl ^x]^ A(1 - log A) where A = -. (8) 

d{q - 1) 

We are now ready to prove Theorem 11.21 Let ^ be a positive real number and define vq = log^ d — ^. 
Applying Formulas (|6]l and (O with x ^ vq + {i + ^) and x — vq — {£ + ^) respectively, we find: 

^ Vo + {£ + ^)] q-' 

and F[Vl ^vo-{e+ i)] ^ • q-' ■ (l - log (^) + ( ■ log q 

^ -^-q-' -{l + i-hgq). 

q - 1 

Theorem 1 1 .21 follows by adding these two inequalities. Corollary [T3] can be now easily deduced. Indeed, 
note that the function v M> IE((Vl — v)"^) is maximal when v is equal to the expected value of Vl and the 
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value taken at this optimal point is the variance of Vl- It is then enough to bound the expected value of 
(Vt — vq)^, which can be done as follows: 



P[{Vl - vo)^ x] -dx 

F[(Fl-wo)' ^ (^+5)'] ■{2£ + l)-de 

^— • / q-^ ■{2 + i- logq) ■ {2£ + l)-de 
9-1 Jo 

q " " 



q-1 \\ogq log^q 



1 

4 
1 
4 
1 
4 



The standard deviation of Vl is then always less than a{q) — yj + ■ ( + q )- The function a 
is decreasing on [2, 00) and then bounded from above by its value at 2 (which is < 6.5). Note furthemore 
that when q goes to infinity, a{q) = ^ + 0{j^). 

1.3 Generalization to block LU decomposition 

Let d = {di, . . . , c?,.) be a tuple of positive integers such that di + ■■■ + dr = d. By definition, a block 
LU decomposition of type d of a matrix M G Md{K) is a factorization M — Ld{M) ■ Ud{M) where 
Ld{M) and Ud{M) are respectively block unit lower triangular and block upper triangular with respect to 
the partition d: 



LdiM) 



fid, 



V 







and UdiM) = 







where the s-th block has size ds and, for an integer n, /„ denotes the identity matrix of size n. Of course, a 
block LU decomposition of type ( 1 , 1 , . . . , 1 ) is nothing but a standard LU decomposition and every matrix 
M S MdiK) admits a block LU decomposition of type (d), whichis simply Af = Id-M. As in the standard 
case, a LU decomposition of type d is unique (when it exists) — which justifies the notations Ld{M) and 
Ud{M) — and, an invertible matrix AI admits such a decomposition if and only if, for all i e {1, . . . , s}, 
its c?i-th principal minor does not vanish. For d as before and oj £ ft, we let VL.d denote the opposite of 
the smallest valuation of an entry of Ld{uj). This defines a random variable VL.d : — N U {00} for 
each d. The aim of this subsection is to study them. Following the same strategy as in the standard case 
(i.e. d— (1, . . . , 1)), our first task is to establish a link between VL,d and the random variables Vi,j defined 
in ilLL2l To shorten notations, we set Is — {di + • • • + ds_i + 1, . . . , di + • • • + dg} and recall that if 
AI g Ald(K) and /, J C {1, . . . , d}, we denote by AIj.j the submatrix of AI consisting of entries whose 
row index and column index are in / and J respectively. For all s € {1, . . . ,t}, we further introduce: 



Corollarv l 1 .5l learns us that the V^.s's are mutually independant for s varying between 1 et r (and d remains 
fixed). The following Lemma shows that their laws are also precisely known. 

Lemma 1.9. For all s G {1, ■••,?'} and all integer v, we have: 

nVd,s ^v] = {l- - q-"-^) ■■■{I- q-"-'^). 

Proof. Throughout this proof, we set a = di + ■ • ■ + dg-i, b = di + . . . + dg and, for i G {1, . . . , ds} 
and i' ~ b + 1 — i, Wi = mm{Vi',i', Vi'^i'+i, . . . , Vi',b). It follows from Corollarv l LSI that Wi follows a 
geometric law of parameter (1 — and furthermore that the Wi's (1 ^ i ^ dg) are mutually independant. 
For all ^ G {1, . . . , ds}, define moreover Se = W\ + • • ■ + Wi. Clearly Sd, — Vd,s- We will prove by 
induction on the couple (£, v) (lexicographically ordered) that: 

P[^, ^ = (1 - - 9""-') ■•■(!- q-''-'). 
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For £ = 1, the statement is true. Assume now that it is true for all {£' , v') with £' < i or i' = i and v' < v. 
The strategy is to decompose the event "Si ^ v" in two parts according to the vanishing or the nonvanishing 
of Wd^^e+i- Clearly, if Wd^-e+i = 0, we have Si = Si-i. On the other hand, if we know for sure that 
Wd^ -i+i does not vanish, one can subtract 1 to it and get this way a new random variable which still follows 
of a geometric law with the same parameter. Hence, one can write: 

V[Si ^v] ^ P[Vt^<i,-^+i = 0] • P[5,_i i^v] + ¥[Wd,^i+i > 0] • ¥[Si i^v-l] 
= (1 - q-') ■ P[5,_i ^ t;] + q-' ■ Y[Si ^v-l]. 

Replacing T\Si-i ^ v] and V[Si ^ t; — 1] by their values (coming from the induction hypothesis), we get 
the desired result. □ 

Remark 1.10. Alternatively, one can notice that Vd,8 follows the same law as the variable "determinant of 
a random matrix of size dg" and then conclude by Abdel-Ghaffar's Theorem. Actually the proof we have 
presented above is very inspired by Abdel-Ghaffar's one. 

Proposition 1.11. We have max(Ki,i; • • • j ^d,r) — Ui?Xdct) ^ VL,d ^ niax(Vd,i, . . . , Vd.r)- 



Proof. We follow the lines of the proof of Proposition ll.6l To avoid confusion, we begin by letting Tj{uj) 
denote the matrix uj computed by Algorithm[T|(run with uj as input) after the j-th iteration of the main loop. 
Pick some s G {1, . . . , r} and set j{s) — di + . . . + ds. We are going to prove the two following statements 
from which the Proposition will follow directly: 

• the determinant of the square dg x dg matrix Tj(^g^{uj)j^j^ has valuation Vd,s(i^); 

• for all < e {s, r}, we have the identity idli^)/*, /a ' Tj^s) ('^)/s./s = 2^j(s)('^)/t,/a- 

The first assertion is easily proved. Indeed, by construction, the submatrix Tj(g){Lu)^i j(^g'jy ^i j(^g)y is 
lower triangular and that its i-th diagonal entry has valuation mm{Vi^i, Vi^i+i, . . . , To prove the 

second assertion, we first remark that, up to replacing w by w + tt^ for a sufficiently large integer N, one 
may assume that ui is invertible. All the matrices ^^(^■j (oj)/, are then also invertible. Consider the matrix 
L e Md{K) whose i-th column is the i-th column of Tj(g-f (uj) where s is the unique such that i G I^. It is 
apparently lower block triangular with respect to the partition d. Furthermore, noting that, for i E Ig, the 
i-th column of Tjf^g-^ (uj) is a linear combination of the first j{s) columns of uj, we see that L^^ ■ oj is upper 
block triangular. Hence, if D is the diagonal block matrix: 



D = 



the factorization uj — {LD ^) • {DL ^u) is the LU decomposition of type d of uj. Therefore Ldiuj) = 
LD~^. Our claim follows directly from this. □ 

From Lemma [L9l we easily derive that g^" ^ ]P[Vd ^ ^ w] ^ • q^"" . Arguing then as in ill.2l one 
can prove analogues of the results we have shown before concerning the random variable Vl'- the expected 
value of VL.d is equal to log^ s + 0(1), its standard deviation is a 0(1) (where the notation 0(1) stands for 
a quantity bounded by a universal constant which can be made explicit) and, actually, we even have a more 
precise (but also more technical) estimation of its law in the spirit of Theorem ll.2l 



2 LU decomposition over a DVR: algorithmic issues 

LU decomposition is a very basic and important tool when we are doing algorithmics involving matrices, 
and especially matrices over a complete DVR. But unfortunately, on some particular inputs, computing it 
may cause important numerical instability; it is the case for instance if the top left entry of the input matrix 
has a very large valuation (compared to the other entries). The first aim of this second section, is to study 
this phemonemon; more precisely, following the ideas of SjT] we will design a new algorithm to compute 
LU decomposition (see Algorithm|2]l and show that the set of unpleasant inputs for which it is numerically 
unstable is very small. 
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In particular, we may expect that if Algorithm |2] is called as a subroutine by an other probabilistic 
algorithm, it will not never generate important losses of precision. In i) il2.2l and l2.3l we will illustrate this 
idea on a particular example: we will propose a probabilistic stable algorithm (based on LU decomposition) 
whose aim is to compute a basis of a coherent over A]^ (where K is the fraction field of a complete DVR) 
from the knowledge of all its stalks. 

We keep the general notations of SjT] Ist i? be a discrete valuation ring whose valuation u/? : i? ^■ 
N U {00} is assumed to be surjective. Let tt be an element of R of valuation 1. Let k (resp. K) denote the 
residue field (resp. the fraction field) of R and set q — Card k. We recall that ti^; extends uniquely to K and 
that, in a slight abuse of notations, we continue to denote by vr this extended map. We recall also that we 
have set = Md{R) and that this space is endowed with its Haar measure. For cj G fi, denote by Wi(a;) 
the valuation of the i-th principal minor of cj and set W = max(iyi , . . . , Wd)- Thanks to Abdel-Ghaffar's 
Theorem (see Q), the law of the W^'s is known: P[Wi w] = (1 - - (J"""^) •••(!- 

for all w ^ and i e {1, . . . , d}. From this, we derive P[VFi > w] ^ |— j- and then: 

V[W >v]^d - ^—^ (9) 

for all nonnegative integer v. Adding all these probabilities, one finds W^W] — log^ d + 0(1) where, as 
usual, the notation 0(1) refers to a quantity bounded by a universal constant. 

2.1 Loss of precision in LU decomposition 

By Formula we know that the entries of L{M) can be all expressed as the quotient of one minor by one 
principal minor. Noting that if x and y are both known with precision 0(7r^) and if y has valuation v, the 
quotient ^ is known with precision at least 0{'k^~'^'"), one may expect that a good algorithm computing 
the LU factorization of lu would shrink the initial precision by a factor 7r2'^(*^) . 

Unfortunately, a quick experiment shows that the naive algorithm based on usual Gauss elimination 
generates losses of precision much more important than that. For example, on a random input matrix 
M e A/25(Z5) given with precision 0(5^), it outputs a matrix L which is in average known up to precision 
0(5^^*^) where c ~ 10 whereas the mean value of 2 ■ W{M) is only ~ 2 • log^ d = A. For matrices of size 
d = 125, the deviation is amplified: we find c ~ 50... to be compared to 2 • log^ d = 6. 

2.1.1 A first simple solution 

Our starting remark is the following: it follows from Cramer like formulae ^ that if M are M' are two 
matrices in Md{R) congruent modulo tt^ (for some positive integer N) such that W{M) < N, that 
W{M') = W{M) and 

L,,,(M) ^ L,,,(Af') (mod n^-^'W)) 

for alH, j G {1, . . . ,d} with i > j. In particular, under the previous assumptions, we have L{M) = L{M') 
(mod 7r^^2W(M)-j jj^j^ j-gsult suggests the following method to compute L{M) with a correct precision 
when M is a matrix known with precision 0{tt^): 

• we lift AI to a matrix M' known with precision 0{'k^ ) for some N' > N; 

• we compute W{M') and L{M') with our favorite algorithm (e.g. Gauss eliminationj^; 

• we answer L{M) = L{M') + 0(7r^-2.w'(M')). 

By what we have said before, our answer L{M) is always correct. Furthemore, if N' is sufficiently large, 
then L(M') will be known with precision at least 0(7r^"2 M'(M')-) ^nd L{M) itself will be known with 
precision 0(7r"-2.w(A/'))^ 

It then remains to find a suitable value for N'. Of course, it will strongly depend on the algorithm we use 
to compute L{M'). Let us study a bit the case of Gauss elimination. Since the successive pivots appearing 
during the elimination have valuations Wi{M'), W2{M'), . . . , Wd{M') and since we are only dividing by 
pivots, the maximal loss of precision is boundedfrom aboveby 2- (Wi(M') + - • ■ + Wd{M'y). In other terms, 

^Generally, these two computations can be done simultaneously. It happens in particular if one uses Gauss elimination. 
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using Gauss elimination, one can certainly compute L{M) with precision 0{ti^ 2 (vKi(m )h \-Wa{M ))y 
As a consequence, it is enough to choose N' so that: 

N' -N^2-{Wi{M) + --- + Wd{M) - W{M)). 

However, at the very beginning, we have not computed the W^j(A/)'s yet. So we cannot figure out at this 
moment what is the best value of A^' (i.e. the smallest one satisfying the above inequality). Nevertheless, 
we know that in average Wi{M') ~ ^ and W{M) ~ Xog^d. To begin with, we can then try to take 

N' = N + \^~\ and see what happens: we do the computation with this particular N', we determine the 
Wi(M)'s, if the above inequality is fulfilled, we are done, otherwise, we determine the right N' and redo 
the computation. Actually, it could happen — but it is very rare — that the first precision 0{tt^ ) does not 
allow us to determine some of the Wi{Mys; in that case, we just guess a new larger N', try with it and 
continue like this until it works. 

Let us finally analyze the complexity of this method in the favorable case where A^' = + ^ is enough. 
In order to fix notations, we assume moreover that doing basic operations (i.e. additions, substractions, 
multiplications and divisions) in R with precision tt^ requires 0((A^ log g)") bit operations where a is 
some constanjl, necessarily greater than or equal to 1. Since the complexity of Gauss elimination is 0{(P) 
operations in the base ring, our method needs: 

0(d3.(iv + |)".log"g) 

bit operations. If d ^ qN, it is quite nice. However, if the opposite situation when d » qN, the dominant 
term in the above complexity is d'^^", which is very large and actually not really acceptable for many 
practical applications. 

2.1.2 A stable algorithm to compute LU decomposition 

In this subsection, we propose and study a different method to compute LU decomposition which has the 
advantage of not requiring to increase the precision at any time and whose complexity is comparable to 
Gauss elimination. Our algorithm is strongly inspired by the constructions of fJT] and especially those of 
^TT2l Here is it: 

Algorithm 2: A stable algorithm to compute the L-part of the LU decomposition 

Input: : A matrix AI of size d x d known with precision 0{7t^) 
Output: : The L-part of the LU decomposition of M 

Notations: : ★ d is the dimension of the matrix M 

: * Aij denotes the {i, j)-th entry of a matrix A 
: -k Aj denotes the j-th row of A 

1 uj^ M; 

2 L ■(r- identity matrix of size d x d; 

3 for j from 1 to d do 



for ifrom 1 fo j — 1 do 

if VR{uji,j) < VR{uji^i) then swap uij and w^; 

iiuji^i ^ then s <— lifted to precision 0{tt^); luj <— luj — s ■ uJi 
for ifrom j + ltodAo Lij ^ ^ + O(7r^-"-'°^'^(0'^«("^>^)-^«("--^))) 



9 return L: 



A first important remark related to Algorithm|2]is the following: at each step, all entries of ut are known 
with precision 0{tt^). Indeed, lu itself is updated only on line |6] and the corresponding computation does 
not affect the precision (because s has been lifted modulo vr^ previously). 

'in usual situations, one can take a = 1 + e for all positive real number e. 
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Correctness of Algorithm |2] We fix an integer j e { 1 , . . . , d} and focus on the matrix lu computed by 
the algorithm after the j-th iteration of the main loop. It is clear that it is obtained from M by performing a 

dot LJ^''^ 

sequence of elementary operations on its j first columns. Thus, for all i > j, we have L{M)ij — j-jf- 

' dct UJ J J 

where / = {1, . . . , j — 1, i} and J = {!,..., ]}. On the other hand, by construction, uj'f\ and Wj are 
two upper triangular matrices modulo tt^ . Their determinants are then congruent to the product of their 
diagonal entries modulo tt^ . Therefore: 

Of course, the value of this quotient is up to some precision. To compute this precision, it is easier to 

work with relative precision (i.e. the difference between the absolute precision and the valuation); indeed, 
we know that the relative precision of a quotient is equal to the minimum between the relative precisions 
of the numerator and the numerator In our case, if we set v = vr{uj[''1) + • • • + and w ~ 

VB.(uj\^j) — vii{u}^''j), the relative precision of the numerator (resp. the denominator) is N — {v + w) (resp. 
TV — v). Thus, the relative precision of the quotient is N — v — max(0, w) and its absolute precision is then 
N — V + min(0, w) (since its valuation is w). The value Li,j computed by Algorithm|2] together with its 
precision, are then correct. 



Precision issues Keeping the previous notations, one certainly have w ^ —v and then iV— w+min(0, w) ^ 
N — 2v. In other words, the (i, j)-th entry of the matrix L returned by the Algorithmic is known with pre- 
cision at least 0(p^~^^'^^^) (recall that Vj{M) denotes the valuation of ujjj at the end of the j-th loop, 
i.e. our previous v). The maximal loss of precision is then bounded above by 2 • max(yi (M), . . . , Vd{M)). 
By the results of ijlj we know that the mean of this upper bound is close to 2 • log^ d, that is the value we 
expected. 



2. 1 .3 Algorithm |2] and Hermite normal form 

Let us denote by H'{M) the matrix lo computed at the end of the execution of Algorithm |2l It worths 
remarking that H'{M) has a lot of things to do with the Hermite normal form of A/. Let us first agree on 
the definition of the Hermite normal form of M: throughout this paper, it will refer to the unique lower 
triangular matrix whose diagonal entries are powers of tt and which is right-equivalent to M (it means that 
H{M) is obtained from A/ by multiplying on a right by a unimodular matrix). We will denote it by H{M). 

Proposition 2.1. Let M G Md{R) known with precision tt^ . We assume that all diagonal entries ofH'{AI) 
are not congruent to modulo and, for all j £ {!,..., d}, we write Hj j{M) = p^^Uj where Vj is a 
nonnegative integer and Uj is a unit. For all i, j £ {!,..., d}, we then have: 

ifKj: H,,j{M) = 

ifi > j : Hij{M) = uj^ ■ Hlj{M) (mod tt^^^O- 

Remark 2.2. Keeping the notations of the Proposition, it is clear that Uj is only known modulo w^^'"^ . The 
congruence of the Proposition is then, by nature, the best one can expect. 

Proof. Let W'{M) be the matrix W' computed by Algorithmic] One can easily check that W'{M) is 
unimodular and moreover that H'{M) — M ■ W'{M). Consequently the Hermite normal form of M is 
equal to the Hermite normal form of H'{M). 

On the other hand, we know that H' [M] has a very particular shape: firstly, it is lower triangular modulo 
TT^ and secondly, by assumption, its diagonal entries are not divisible by tt^ . Thus, H{M) is obtained from 
H'{M) by clearing one by one its entries lying above the diagonal and by dividing its j-th column by Uj. 
But, if Hij{M) = TT^Vi,j (for some pair with i < j), one clears the (i, j)-th entry of H'{M) by 
doing the following elementary operation on columns: Hj{M) -s— Hj{M) — TT^-^u-^v.^jHliM). Hence 
clearings do not affect the value of H- j{M) modulo Tr^^"^. The Proposition follows easily from this 
observation. □ 



11 



2.1.4 The notion of L'V decomposition 



The L-part of the LU decomposition has of course very nice abstract properties but unfortunately does not 
behave very well regarding to precision. Indeed, as we have seen before, if a matrix M is known modulo 
TT^, it is not true that L{M) is known with the same precision. But, beyond that, the precision data attached 
to L{M) is not uniform in the sense that all entries of L{M) are not known with the same precision. In 
order to tackle this problem, we introduce the following definition. 

Definition 2.3. Let M G Md{R) and be a positive integer. A L'V decomposition of M modulo TT^ is 
a couple of d X d matrices {L', V') such that L' = MV' (mod tt^) and L' and V' are lower triangular 
modulo TT^ and upper triangular modulo tt^ respectively. 

If there exists a diagonal entry of L' which is congruent to modulo tt^, (L', V') is said to be degener- 
ate. Otherwise, it is nondegenerate. 

Remark 2.4. It is easy to see that if {L',V') is nondegenerate, then all diagonal entries of V are not 
congruent to modulo tt^ as well. 

It is not difficult to modify Algorithm |2] so that it computes a L'V decomposition modulo tt^; we end 
up this way with Algorithmic] On the other hand, it is worth noting that L'V decomposition is closely 



Algoritlim 3: An algorithm to compute a L'V decomposition 



Input: : A matrix M of size d x d known with precision 0{tt^) 
Output: : A L' V decomposition of M modulo tt^ 

1 w ^ M; 

2 L',V' ■(r- two new matrices of size d x d; 

3 W' ^ identity matrix of size d x d; 

4 for j from 1 fo d do 
for ifrom 1 fo j — 1 do 
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if (wi.j) < VB.{i^i,i) tlien swap Uj and uji, swap W'j and W^; 
ituji^i ^ tlien 



lifted to precision 0{'k^)\ 
- Lo, ~ s ■ uj,: W'^W'-s- W': 



11 return L', V; 



oj f, V\ ^ W'f, 



related to LU decomposition. The following proposition makes this statement precise. 

Proposition 2.5. Let M e Md{R), N be a positive integer and {L' , V') be a nondegenerate L'V decom- 
position of M modulo TT^. Then M admits a LU decomposition and for all (i, j) with 1 ^ i < j ^ d, one 
have: 

L,j{M) = (mod T,N-v,-m^^(0,vnmj-vn{Ll^))-^ 

with Vj = ^RiK.k) - VRiKk)- 

Proof. Left to the reader (the arguments are very similar to those detailed in i j2.1.2l i. □ 

Of course, executing first Algorithm[3]and then applying the result of Proposition l2.5l is almost the same 
than running directly Algorithmic Nevertheless splitting Algorithm |2] in two parts can be very useful for 
some applications (we will see an example of this in i j2.1.5b because, as we have already said before, the 
pair {L' , V') is generally easier to manipulate than L{M) since it carries a flat precision (and, in addition, 
it consists of two integral matrices if M is itself integral). 

2.1.5 Complexity and Hafner-McCauley's algoritlim 

It is easily seen that the asymptotic complexity of Algorithmic is 0{d^) (operations in the base ring R) 
where d denotes the size of the input matrix. It is then similar to the complexity of usual Gauss elimination 
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whereas it is true that our Algorithm|2]runs a Uttle bit more slowly because it basically makes more swaps 
and copies. 

When precision is not an issue (e.g. when we are working over an exact ring), Hafner and McCauley 
showed in |5j how to reduce the computation of the LU decomposition to matrix multiplication and got 
this way a nice recursive algorithm that computes the LU decomposition of a matrix in only 0((i") oper- 
ations where u is the exponent for matrix multiplicatioifl The aim of this subsection is to extend Hafner- 
McCauley's algorithm in our setting where we want to take care of precision. 

A preliminary result about Algorithm |3] Roughly speaking. Algorithm [3] clears the entries of uj lying 
above the diagonal in the colexigographic order. We would like to study what happens if we decide to clear 
these entries in a different order 

Definition 2.6. Let a < 5 be two positive integers and set /a,b = {(i, j) G | a ^ i < j ^ b}. A total 
order =:<; on la^b is nice if: 

• for l^i^i'<j^d, one always have =^ (i' and 

• for all and G la.b such that j ^ i', one have ^ {i' 

Remark 2.7. It is easy to check that the colexicographic order on la.b is nice. However, it is not the only 
one: the lexicographic order, for instance, is nice as well. One can also build recursively nice orders on la.b 
as follows. Fix an integer c between a and h and pick and =<;2 two nice orders defined on la.c and Ic+i,b 
respectively. Consider also a third order =^3 defined on the cartesian product {a, . . . , c} x {c + 1, . . . , fc} 
and satisfying the first condition of Definition 12.61 Now define a new order ==; on la.b by agreeing that 
Ia,c ^ {a., • ■ • 7 c} X {c + 1, . . . , &} =^ ^c+i.Jl and furthermore that =^ agrees with =^1, =^2 and =^3 on la^c, 
Ic+i,b and {a, . . . , c} x {c + 1, . . . , 6} respectively. A quick check then shows that =^ is nice as well. 

If ^ is a nice order on Ii.d = {{hj) G | 1 ^ i < j ^ d}, let us agree to use the expression "to 
execute Algorithm |3] with respect to to mean that we execute this algorithm but, instead of running 
through all (i, j) e Ii.d according to the colexicographic order, we run through these pairs according to ^ 
and execute line [TO] when i = j — 1. 

Proposition 2.8. When they are called on the same input, Algorithm\3\and Algorithm\3\executed with 
respect to a nice order return the same answer 

Proof. Easy check. □ 

Description of the algorithm Suppose that we are given a matrix M £ Md{R) known with precision 
0(7r^). The basic idea (which comes from Hafner and McCauley) is to obtain a recursive algorithm to 
compute the LU decomposition and, doing this, to replace as much as possible elementary operations on 
rows by matrix multiplication. Moreover, in order to avoid many problems related to precision, it would be 
really better to work with L' V decomposition instead of LU decomposition. Actually, for the purpose of 
the recursion, we will not just need the matrices L' and V' but also H' (which is the matrix uj at the end of 
the execution; see i )2.L3l l and W' . The prototype of the algorithm we want to design is then: 

LV : Af h-^ {L',V',H',W'). 

Proposition l2.8l together with the recursive construction of a nice order detailed in Remark l277l su g ge sts the 
following strategy for a recursive implementation of LV: 

L we start the computation of a L'V decomposition of M but stop it after d' columns for some d' < d 
(e.g. d' = [f]); 

2. we clear all the entries in the d' x (d — d') top right corner of the matrix uj we have ended up after the 
first step; 

■^Nowadays, the best known value for uj is 2.376 but, unfortunalety, the corresponding algorithm due to Coppersmith and Winograd 
(see [3|) is not efficient in practice (even for very large d) because the constant hidden in the O is quite large. A good compromise is 
to use classical Strassen's algorithm whose asymptotic complexity is a little bit worse — exactly 0(d'°S2 7) _ but which is easy to 
implement and works very well in practice. 

""By this inequality, we mean that elements in la.c are all less than those in {a, . . . , c} X {c-\- 1, . . . ,b} and, in the same way, that 
the latter elements are less than any pair in I^+i j,. 
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3. we finally compute a L' V decomposition of the {d — d') x (d — d') bottom right corner of lo. 



It turns out that the first step can be computed in a recursive way. Precisely, we decompose M as a block 
matrix 



Af = f \ (where Mi has size d' x d' 

\ M^ M4 J 



we call recursively the routine LV on the input M\ and then recover the matrices L', V and uj (as they have 
to be just after the first step) using the following formulas: 



M3 • Wi Mi 



where the quadruple [L'l^ V(, H'^, W{) is the output of the recursive call of LV. Last but not least: remark 
furthermore that, proceeding this way, we are replacing elementary operators (on the columns of by 
matrix multiplication (by V{ and W{). It is exactly the benefit we were looking for! 

Let us now focus on step 2. With the notations above, it consists in clearing aU the entries of M2 
(using eventually the diagonal entries of H'l as pivots). Of course, this can be done just by running the 
corresponding part of Algorithm[3] Nevertheless, we do not want to proceed exactly along these lines but 
we would like instead to use a recursive version of this algorithm in order to take advantage again of the 
complexity of the matrix multiplication. Writing such a recursive version is actually very similar to what 
we have done before. In order to have more coherent notations, let us rename H'l and AI2 to X and Y 
respectively and write: 

■Xi \ /Yi Y2 



--it I) ^ -(S I)- 



Note that Xi and X/i are then lower triangular modulo tt^. We can then proceed recursively along the 
following lines: 



1 . we clear Yi using Xi as pivot; 

2. we clear Y2 using the new Xi as pivot; 

3. we clear ¥3 using X4 as pivot; 

4. we clear 14 using the new X^ as pivot. 



Each of these steps can be done recursively. As in the previous case, we just need to be careful and let 
our recursive routine return not only the new matrix X gotten after clearing Y but also the transformation 
matrix T such that {X Y ) ■ T = {X' ) where X' is the new X mentionned previously. Indeed, this 
matrix is needed to update X and Y after each step. 



A brief study of complexity Let us denote T'{d) the complexity of the clearing algorithm we have just 
described (i.e. the number of elementary operations on R performed by this algorithm when the size of 
the input matrices is d) and by T{d) the complexity of our complete recursive algorithm computing a LV 
decomposition. From the description of these algorithms, we find: 

T[d) = 2-T(f)+r'(f) + 0(d") (10) 
T'{d) = 4-T'(f) + 0(d") (11) 

where we recall that lo is the exponant of the complexity of matrix multiplication. Since a d x d matrix 
have d^ entries, one certainly have lu ^ 2. For simplicity, we assume that ui > 2 (we recall that the 
fastest asymptotic algorithm known today corresponds to w ~ 2.376). It is then a classical exercise to 
deduce from the recursion formula (fTTT i that T'{d) = 0{d'^). Knowing this, equation ( fTOb becomes T{d) = 
2 ■ r(f) + 0(d'^) and then yields T[d) = ©(d'^) as expected. 



14 



Algorithm 4: Computing the L-part of the LU decomposition of type d 
Input: : A partition {di, ... ,ds) of a positive integer d 

: A matrix M G Md{R) known with precision 0{'k^) 
Output: : The matrix Ld{M) 

1 w M; 

2 L <~ zero matrix of size d x d; 

3 5^1; joi~ 0; 

4 for j from 1 to d do 



5 


for i from 1 fo j — 1 do 


6 




if vji{ujij) < vji{uji^i) then swap ujj and cj^; 


7 




if Wi^i 7^ then s hfted to precision 0(7r^); luj <— loj — s ■ w^; 


8 


ifi = jo + ds then 


9 




for j'/'w'w jo + lto jo + ds do ^ —J— • ; 
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for j'from jo + 1 to jo + 4 do Lj, ^ Lj, - E?=/+i j' • -^^'i 
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12 




for j' from jo + 1 to jo + ami i' from jo + 1 to d do set precision of Li' j' to 0{ 
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jo ^ jo + ds; + 



14 return L; 



2.1.6 Block LU decomposition 

The results of i l2. 1.21 extend to block LU decomposition using i ll. 31 Indeed, a close look of the proof of 
Proposition 11.111 shows that one can compute the block LU decomposition of type d = {di, dr) of a 
matrix M e Md{R) using a slight modification of Algorithmic which consists in updating the matrix L (on 
line[8]l only if j is equal to some di + ■ ■ ■ + ds and clearing the entries of L below the diagonal of the s-th 
block just after this update (cf AlgorithmH]). Furthermore, if the input M is known up to precision 0{'k^), 
the precision of the matrix L returned by Algorithm 2] is at least 0(7r^~^'^^ i^*^*^'). In average, the loss of 
precision is then bounded by 2 ■ £[1/^.^] ~ 2 • log^ s. 

All other results proved previously for classical LU decomposition (relation with Hermite normal form, 
notion of L' V decomposition, Hafner-McCauley's improvement) also extend almost verbatim to block LU 
decomposition. We will not explain it in details here (but let the exercise to the reader). 

2.2 Simultaneous PLU decompositions 

As we have already said before, a LU decomposition may fail to exist for some particular matrices. Nev- 
ertheless, it is well known that all matrices over a DVR admit a PLU decomposition (recall that a PLU 
decomposition of a matrix M is a factorization M = PLU where P is a permutation matrix and L and U 
are as before) and, in general, that several matrices P are possible. 

Assume now that we do not pick just one matrix, but a (finite) family of matrices (Mi, . . . , Af„). The 
question we would like to address is the following: does there exist a "simultaneous PLU decomposition" 
of the A/i's, that is PLU decomposition of each Mi with the same matrix P. If we want as before P to be a 
permutation matrix, the answer is negative in general. However, if we relax this condition and require only 
that P is invertible (which is enough for certain applications, see i l2.3l for a concrete example), the answer 
is positive (at least if the base field is infinite). 

The aim of this section is to study this notion of "simultaneous (block) PLU decomposition" over a base 
field which is the fraction field of a discrete valuation ring. 

2.2.1 The basic result 

Let {Ml , . ■ . , Mn) be a family of square d x d matrices over K and fix a partition d — {di, ... ,ds) of d. 
There exists an obvious probabilistic algorithm to compute a simultaneous block PLU decomposition (of 
type d) of the M^'s: we choose P at random and compute the block LU decomposition of the P^^Af^'s. 
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The aim of this paragraph is to prove that this algorithm works quite well in the following sense: not only 
it finds very quickly a matrix P that does the job, but it furthermore finds quickly a matrix P for which all 
entries of P, P~^ and the i^'s are kwown with a good precision and do not have a too small valuation. Our 
precise result can be stated as follows. 

Theorem 2.9. Let n be some positive integer. Suppose that for all m € {I, ... ,n} we are given a matrix 
Mm G Md{K) together with a finite sequence = (dm.i, • ■ • idm,r,n) of positive integers whose sum 
equals d. Let e be a real number in (0, 1) and take v an integer ^ log^ C'^"^'^'!^'"'"" ) — log^ e. Then, a random 
matrix a; G satisfies the following conditions with probability at least 1 — e: 

• u! is invertible in Md{K) andu!~^ G Md{R); 

• forallm€ {1, ... ,n}, the matrix ujMm admit ablockLU decompositionof type d^and Ld {ujMm) G 
7r--Md{R); 

Moreover if the M^'s all lie in Md{R), are invertible in this ring and are only known with precision 0{tt^), 
one can furthermore require (without changing the probability) that, on each input ujMm, Algorithm^ 
outputs Ld^ {ujM„i) with precision at least 0{tt^^'^'"). 

Remark 2.10. It is also possible to bound the loss of precision if we drop the hypothesis of inversibility of 
the Mm's. The valuations of their determinants then enter into the scene. (The exercise is left to the reader) 

Proof. Let us begin by proving the first assertion. Replacing eventually Mm by A/m + 7r^ for a large integer 
N, one may assume that all A/^'s are invertible in Md{K). Furthermore, since multiplying Mm (1 ^ m ^ 
n) on the right by an upper triangular matrix with coefficients in R does not change the matrix Ld^ [ujMm) 
for any a; G Md{K), one can even safely assume that all M^'s are invertible in Md{R) = ^. For all 
TO G {1, . . . ,n} and all s e {1, . . . ,r„i}, define W±s,m : -> N U {oo} by Wd^^sM^) = V^,s(wM„) 
(where Vd,s is the random variable defined in m.3\ and set: 

W = ma,x{v Bidet), Wd^s,l, Wd^s,ri-l,Wd^^2,l, • ■ • , Wd^^2,r2-1,- ■ ■ , W^d„,n,l, • • ■ , W^d„,n,r„-l)- 

Since Af,„ is invertible in 57, Lemma [L9l implies that 



nwd„^,m^s ^v] = {i- q-'"-') ■■■{!- g-^-^^o < 1 - 




for all TO, s and v. Furthermore by Adbel-Khaffar's Theorem (see Theorem 1 of 1 1 1) we also know the law 
of the random variable ■yfl(det); we have P[vR{det) w] = (1 + + g^^^^) • (1 + q^''^'^). Let 

us simplify this Formula and just remember that P[ti/{(det) ^ v] ^ 1 — (^^"^^ + g^"^^ + • ■ • + q-v-d) ^ 
1 — l^j. We can now estimate the law of W as follows: 

P[W >v]i^ PK(det) >v]+ J2 ^[^™.^ > ^] ^ ■ ( '"'^"'^j''"'" + ^) < £• 

Proposition 11.111 shows that Ld^^^ (ujMm) e n^^^'^'' Md{R) for all w G 17 and all to G {1, ... , n} and, on 
the other hand, it is clear that uj^^ G tt"*^*^"^ Afd(i?) because vii{det w) ^ W{uj). It is enough to conclude 
the proof. 

The second assertion (concerning precision) is now clear. □ 

Remark 2. 11. One may wonder if the bound log^( ''^'^'^'j''''" ) — log^ e is sharp. Actually, it cannot be for any 
data of (Ml, . . . , A/„). Indeed, an integer v satisfies the required conditions of Theorem l2.9l for the families 
(M, . . . , M) and {d, . . . ,d) (n times) if and only if it satisfies the same conditions for the family reduced 
to the unique matrix M. So if Mi = • • • = M„ and di = . . . = d„, one can certainly improve the bound 
logg(^~r) ~ logg £ by removing the facteur n in the first log. Nevertheless, by using similar methods as 

those of 2] one can prove, first, that the result of Theorem |2.9| fails if v <C logg C"'^^^^^'^^"'''"'' ) — log^ e and, 

second, that if Mi, . . . , M„ are themselves chosen randomly, it even fails for v <^ log^ C'^^'^^''" ) — log^ e 
(i.e. the given bound is sharp). 
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2.2.2 A slight generalization 

For the application we have in mind (see H2.3\ . we will need a slight generalization of Theorem |Z9] where 
the matrices Mi , . . . , Af„ on the one hand and the matrix lu on the other hand are not defined over the same 
field. Let K he a finite extension of K. A classical result asserts that the valuation vji extends uniquely 
to K. Let R be the ring of integers of K, that is the subset of K consisting of elements with nonnegative 
valuation. Set finally: 

oo 
i=l 

It is easy to check that q ~ 1 — < Il{q) < q — 1. 

Theorem 2.12. Let n be some positive integer. Suppose that for all m G [1, . . . ,n} we are given a matrix 
Mm G Md{K) together with a finite sequence — (rfm.i, • ■ • idm.rm) of positive integers whose sum 
equals d. Let e be a real number in (0,1) and take v an integer ^ log^ ( ) — log,; £■ Then, a random 
matrix a; G f2 satisfies the following conditions with probability at least 1 — s: 

• u! is invertible in Md{K) andoj^^ G tt^"" Md{R); 

• forallm£ {1, ... ,n}, the matrix ujMm admit ablocJcLU decompositionof type dand Ld [uiMm] £ 
TT-^MdiR); 

Moreover if the Mm's all lie in Md{R), are invertible in this ring and are only known with precision 0{tt^), 
one can furthermore require (without changing the probability) that, on each input ujMm, Algorithm^ 
outputs Ld^ {ujMjn) with precision at least 0{tt^~'^'"). 

Remark 2.13. Once again (see Remark |2.10| |. one can bound the loss of precision as well if we drop the 
hypothesis of inversibility of the A/„j's and put into the machine the valuations of all det Af,„. 

We now start the proof of Theorem 12.121 it will occupy the rest of this subsection. As in the proof 
of Theorem 12.91 we start with the first assertion and assume that Mm is invertible in Md{R) for all m. 
However, in our new settings, this fact no longer implies that ujMm runs over Md{R) when uj runs over 
Md{R). Thus, we can no longer work with the random variables lu ^ Vd ^iiojMm) and we need to 
modify a bit our strategy. Actually, since we just want to bound from above — and not from below — the 
valuation of the matrices Ld^^ (ujAIm), we can argue first assuming that m is fixed and then add probabilities. 
Moreover, by the proof of Proposition 11.111 (see also Formula ^ when d,„ = (1, . . . , 1)), bounding the 
valuation of Ld^ (wMm) reduces to bounding the valuation of the {dm,i + . . . + (im.s)-th minor of LuAIm 
for all s G {1, . . . , — 1}. Thus we first fix to G {!,...,«} and s G {1, . . . , — 1} and look for 
an upper bound for the valuation of the jm{s)-iti minor of ujMm where, by definition, jm{s) = dm,i + 
■ • • + dm.s- For 1 ^ i ^ jm{s), we are going to define a random variable Wi : ^ U {oo} 
where e is the ramification index of K/K (i.e. the index of Vfj{K*) is vu{K*)). The construction of the 
Wi's is achieved by applying the classical algorithm of LU decomposition. We pick w G and first set 
M^^^ = (wMm){i....,r,„}.{i,....r„j G Mr{R). Let j be the first indcx for which '!;fl;(M.[^j' ) is minimum 
among the valuations of all entries of the first line of M'^^\ Let Af*^^-' be the matrix obtained from A/^^^ by 
swapping the j-Xh column with the first one and by clearing all the entries of the first row (expect the first 
one) by pivoting, i.e. adding to each column (expect the first one) a suitable multiple of the first one. The 
matrix 

Af (2) looks like 

/★ •■• 0\ 

* 

V */ 

We now continue this process: we select the first index j ^ 2 for which vr{M2 ■) is minimal, we obtain 
Af by putting the column in the second position and clearing all the other entries on the second row. 
Repeating this again and again, we obtain a finite sequence Af . . . , Af '^'^ of matrices and the last one 
is lower triangular For i G {1, . . . , jm(s)}, we define Wi(a;) as the valuation of the i-th diagonal entry of 
M(-''"('')) (or equivalently of M'^^'^ for some i ^ i). It is clear that the determinant of the j„i(s)-th principal 
minor of ooMm is equal to Wi{uj) + W2{lo) + • • • + Wj^^i^g){uj). We need to determine the law and the 
correlations between the Wi's. We begin by a Lemma. 
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Lemma 2.14. Let J : B!^ ^ FT be a surjective map. Then 
for all nonnegative integer v. 

Proof. Let k denote the residue field of R; it is a finite extension of k. Since / is surjective, it induces a 
surjective A;-linear map f : k'^ ^ k'^ over the residue field. Moreover, the image of / is generated over k by 
f{k'^)- Thus dinifc /(fc'^) = dim^. /(fc'') = r. This fact implies the existence of a fc-linear map g : k^ k^ 
such that the composite g o f : k'^ k^ is surjective. Let g : ^ R^ be any i?-linear lifting of g. 
The i?-linear morphism h = g o f : R'^ ^ R^' induces a surjection over the residue field and thus is itself 
surjective. Furthermore, it is clear that h{x) is divisible by tt" if f{x) is. Hence, for x staying in R'^, we 
have P[/(a;) S tt'"R^] ^ P[h{x) € 7t'"R^] and we are reduced to prove the Lemma with / replaced by h. 
(In other words, we may assume that K — K.) 

By the structure Theorem for finitely generated modules over a principal domain (recall that i? is a 
principal domain), there exists a basis (ei, . . . , Cd) of i?'' such that the first [d — r) vectors ei, . . . , Cd-r 
form a basis of ker/i. Now, using that h is surjective, we easily see that a vector x = X]i=o ^ -^'^ 
satisfies h{x) G 'k'"R^ if and only if Xi is divisible by tt"" for alH > d — r. But, the probability that such an 
event occurs is g^*"" and we are done. □ 



Corollary 2.15. For all integers wi, . . . , ^^^(s). have: 



— {r-\-l — i)vi 



1=1 

Proof. For uj E and i E {1, ■ • • ,jm{s)}, we denote by uii the i-th row of ut (and consider it as a vector 
of R'^) and by i\f ('^(a;) the matrix defined above. Let F^^^ be the submatrix of consisting of its first 
columns and let /^^^ : R'^ ^J™!") the i?-Hnear map whose matrix is 'F^^^. The fact that Mm 
is invertible implies that /*^^' is surjective. Lemma |2 . 1 41 applied to /*^^' yields: 

V[Wi ^ vi] q-^"^'^''^''K (12) 

Now remember that M*^^^ (uj) is obtained from M*^^' (uj) by performing a sequence of elementary operations 
on columns. It then exists a matrix p(i)(w) such that M'^'^\uj) = M^'^\uj) ■ p(i)(w). Clearly P'^'^\lj) 
depends only on uji and we will denote it P^^H^i) in the sequel. Set F'^'^^uji) = p(^) • P^^^^i) and let 
: R'^ PJ-^") denote the map whose matrix is 'P^^^(a;i). It is surjective and one can then apply 
Lemma l2.14l to the composite pij^^^(s)-i ° /^^H^i) where pr^^(^-)_j^ : pj"'('') — > is the projection 

on the first coordinates. It gives P[W2 ^ V2 \ uji ~ xi] ^ q-Umis)-i)v2 Jqj. ^jj rj.^ ^ j^d integrating now 
over xi and using ( fT2b . we get 

¥[Wi ^ vi and W2 ^ V2] ^ g-J-^")"! • 

The Corollary follows by repeating jm{s) times the previous argument. □ 

If we denote by Sg.jni^) the valuation of the determinant of the jm(s)-th principal minor of uiMm, 
Corollarv l2. ISI allows us to do the following computation: 

VlA ^"3m(s)='" + l '"iH ^^jm(s)=^ + ^ 

^ q^'^'^l ^ q-(»'2+2'»3+---+0™(s)-l)«j„(,)) 

V2A ^"i„,(s)^"+l 
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It is time now to free s and m: summing the above estimation over all possible s and to, we find that 
= maxs.m (5s,m(w) is greater than v — which implies that Ldi^Mm) does not lie in tt^"" Md{R') — 
with probability at most (ri + ■ • • + r„ — n) • q^^ ■ n(g)^^ and consequently that uj does not satisfy the 
conditions of the first statement of Theorem l2. 121 with probability at most; 

_„ / 1 , ri H h r„ - n \ _„ ri H h r„ 

^ q 



\q-l n(g) J Uiq) 

Hence if v is chosen ^ log^{ ^^^'^^^'^ ) — \og^{e), this probability is less than e: the first part of Theorem 
I2.12l is proved. 

The second part now follows easily: indeed, we know that, on the input ojMm, the Algorithm|2]decreases 
the precision by a factor that cannot exceed tt^ (^m,,^ ) a fortiori, by a factor that cannot exceed 

^25{M) \yjjere 6 is the random variable defined above. The conclusion follows from this. 

Remark 2.16. The bound of Theorem l2.12l is sharp if Mi , . . . , M„ are chosen randomly among all square 
d X c? matrices with coefficients in R. However, it is not true in general and it is even not true if Mi , . . . , M„ 
are chosen randomly among all matrices over R. Indeed, in that case, using results of fJT] one can prove 
that, in average, the better possible bound for v is given by: 

[k : K]-^ ■ (log,(n±j^) - log^(£)) + 0(1) 
with an extra factor [K : K]^^, which can be very small. 



2.3 Modules over K[X] and sheaves over A\ 

Let X denote an affine curve over K and A — K[X] be the ring of regular functions over X. It is well 
known that the category of coherent sheaves over X is equivalent to that of finitely generated modules over 
A. In particular, the data of a submodule M C A"^ (for some fixed integer d) is equivalent to the data of a 
coherent subsheaf AA C Oj^. Nevertheless, these two objects are of different nature and we would like to 
represent them in two different ways: 

• a submodule M C A'* by a matrix of generators 

• a subsheaf Ai C Oj^ by the data of the stalk A4x C ^ for each closed point x G X (note that 
this inclusion is not trivial for only a finite number of points x). 

Since these objects are supposed to be equivalent, it is natural to ask if one can find an efficient way to go 
from one representation to the other Actually going from the global description to the local one is quite 
easy: it suffices to localize at each point x. Contrariwise, going in the opposite direction is not so obvious 
and will be discuss now. 

From now on, we assume for simplicity that X is the affine line (and leave to the reader the exercise 
to extend our constructions to a more general setting). With this extra assumption, the ring A is nothing but 
the ring of univariate polynomials with coefficients in K. 



2.3.1 Rephrasing our problem in concrete terms 

For all irreducible polynomials P e if [^], let Ap denote the completion of A for the P-adic topology, that 
is Ap — ^im^ A/P"^ A. Concretely Ap can be identified with a ring of power series with coefficients in the 
residue field Kp — A/ PA in one indeterminate Xp. This variable Xp should be thought as "X — ap" 
where ap is a (fixed) root of P in Kp. Under the identification Ap ~ iirp[[Xp]], the natural embedding 
A Ap is just the Taylor expansion at ap: 

1=0 

Let Pi, . . . , P„ be the minimal polynomials of oi . . . , a„ respectively and, for simplicity, set Am = Ap^. 
The question we have addressed earlier is then equivalent to the following: given, for all to G {1, ... , n}, a 
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submodule Mm ^ti ^ee of maximal rank, how can one find expficitely a A-module M C A'^ such that 
Am <SiM = Mm (as a submodule of A^) for all m and A p <Si M = Aj, for sil other P? 

One can actually rephrase again this question in very concrete terms by taking basis everywhere. Indeed, 
if iJ is A or one of the ^Im's, any free submodule of B"^ of rank d can certainly be represented by a square 
d X d matrix with coefficients in B: the module is recovered from the matrix by taking the span of its 
column vectors. Note furthermore that two matrices G and H defines the same module if and only if there 
exists an invertible matrix P with coefficients in B such that G = HP; if this property holds, we shall 
say that G and H are right-equivalent. Since all our base rings are principal domains, we know that any 
matrix G G Md{B) admits a factorization G = MDN where M and N are two invertible matrices, D is 
diagonal and each diagonal entry of D divides the next one. Up to replacing G by a right-equivalent matrix, 
one can furthermore assume that is the identity matrix, i.e. that G has the particular form G = MD. 
Moreover, if B is one of the A^'s, it is safe to assume that the diagonal entries of D are all some powers of 
the variable Xm since all nonvanishing element of Am can be written as a product of an invertible element 
with a power of Xm- In that case the data of D is then reduced to that of a nondecreasing sequence of 
integers ni < . . . < n^. 

With aU these remarks, our question becomes: 

Problem 2.17. Given for all to, an invertible matrix Mm G Md(Am) and a nondecreasing sequence of d 
integers e^,! ^ • • • ^ em,d, how can one construct explicitely a couple (M, D) of matrices over A such 
that: 

i) the matrix M is invertible in Mci{A) ; 

ii) the matrix D is diagonal and each of its diagonal entry divides the next one ; 

iii) for all m G {1, • • • jfi}, the matrix MD is right-equivalent to MmDm over Am where Dm = 
Diag(X^'-\...,X^'-<^); 

iv) for all irreducible polynomial P € K[X] which is not one of the P^'s, the matrix MD is right- 
equivalent to the identity matrix over Ap. 

2.3.2 The answer 

We consider, for all to G {!,..., n}, an invertible matrix Mm G Mii{Am) together with a nondecreasing se- 
quence of d integers e^.i ^ • • • ^ e-m,d- Our aim is to construct a couple (M, D) satisfying the Conditions 
i), ii), iii) and iv) above. Firstable, we define the matrix D as follows: 

D = Diag(Pr'' • • • P„^-^ , . . . , Pi'-' ■ ■ ■ P>'). 
It clearly satisfies Condition ii). 

Lemma 2.18, Let M G Md{A). 

a) Assume that, for all m G {1, . . . , n}, the matrix M considered as an element of M(i{Am) fvia the 
natural embedding A — ^ Am) is congruent to Mm modulo X^'"^^ . Then, the couple (M, D) satisfies 
Condition iii). 

b) Assume moreover that M is invertible in Mi{A). Then the couple (M, D) satisfies Conditions i), ii), 
iii) and iv). 

Proof. Note that in the ring Am, the polynomial Pm is equal to the product of Xm by a unit whereas all 
other Pto''s (for to' ^ to) are invertible. We deduce from this that D is right-equivalent to Dm over Mm- 
Hence, our first hypothesis impUes that MD is right-equivalent to a matrix congruent to MmDm modulo 
j^e^,d+i other words, there exists a matrix Q G GLd{Mm) such that MD is right equivalent to: 

MmDm + X^-'+^Q = MmDm ■ [h + " Diag(X^-<^-^-^ , . . . , x^'''''-'^'-' , 1) ■ Q] 

where Id is of course the identity matrix. The last factor (the one between brackets) is a matrix over Am 
congruent to identity modulo Xm. It is therefore invertible. It follows that MD is right-equivalent to 
MmDm, and part a) of the Lemma is proved. 



20 



We assume now that M is invertible. Then, clearly. Condition i) holds. Moreover, we have already seen 
that Conditions ii) are iii) are fulfilled. It is then enough to prove Condition iv). Let P G K[X] be an 
irreducible polynomial different from all the -P„i's. All f,„'s are then invertible in Ap and, consequently, 
so is the matrix D. Since M is itself invertible, the product MD belongs to Gl^d{Ap) and is then right- 
equivalent to the identity matrix. □ 

It is actually not difficult to produce a matrix M satisfying the assumption of the Lemma l2.18l a). Indeed, 
the identification ^m/Xm" '* Am — A/ Pm"'''''^^^ shows that the congruence M = Mm (mod X^''^'^^) 
is equivalent to M = M'^ (mod Pm" ''^ ) for a certain matrix S Md{A). Hence, finding a convenient 
M is just a direct application of the Chinese Remainder Theorem (recall that all Pm's are irreducible and 
pairwise distinct polynomials). 

Producing a matrix M satisfying also the second assumption of Lemma l2.18l is a bit more tricky but 
can be achieved using block LU decomposition. For m G {!,..., n}, let r,,, be the numbers of differents 
values taken by the sequence {em,l^ ■ ■ ■ , em,d) and let d„i,s (1 =^ s ^ r„i) denote the number of times this 
sequence takes its i-th smallest value. We then have: 

Assume now for a moment that all M^'s admit a block LU factorization Mm — LmUm of type d,„ — 
{dm,i, ■ ■ ■ , dm,r,^)- Sincc Mm is invertible, so is Um- Let Vm be the matrix obtained from Um by mul- 
tiplying its («, j)-th entry by X^'' ^"^ ' (note that the exponent is always nonnegative when the (i, j)-th 
entry of Um does not vanish). Obviously, Vm. is again upper triangular and its diagonal entries are equal to 
those of Um- Thus [/,„ and Vm share the same determinant and Vm is invertible. Moreover, we check that 
MmDm = LmDmVm, from what we derive that MmDm is right-equivalent to LmDm- Since all L„j's are 
unit lower triangular, there certainly exists a unit lower triangular matrix L G Md{A) which is congruent to 
Lm modulo Xj-i^'''^^ for all m. Such a matrix is apparently invertible and also satisfies the assumption in 
part a) of Lemma l2.18l We can then simply take M — L. 

Now let us go back to the general case where some Mm might not have a block LU decomposition of 
type dm- In that case, we denote by Mm{0) € Md{Km) the image of Mm under the canonical projection 
Am Am/PmAm = Km (or, cquivalcntly. Am — Km[[Xm\] Km)- The coefficients of A/,„(0) then 
aU lie in Km, which is a finite extension of K. We can therefore apply Theorem 12. 121 which implies in 
particular the existence of a matrix lo G GLd{K) such that uj ■ Mm{Q) has a block LU decomposition of 
type for all to. Lemma [2.191 below shows that this decomposition lifts to a LU decomposition of type 

dm ofujMm- 

Lenuna 2.19. Let L be a finite extension of K. Pick M G A/d(L[[y]]) and denote by M(0) its image in 
Md{L) under the projection P[[l^]] — > L, K 0. Assume that M(0) is invertible and admits a block LU 
decomposition of type dfor a certain partition d of d. Then M also does. 

Proof. Write d— {di, . . . , d^) and set as usual j{s) — di + ■ ■ ■ -\- dg for all s G {1, . . . , r}. It is enough to 
check that, for all s, the j'(s)-th principal minor of Af, say 6s{M), is invertible in But, being 

a local ring, 5s [M] is invertible if and only if its image Sg (Af ) (0) is. Now remark that this image is nothing 
but the corresponding minor of M{{)): in other words (5j(M)(0) ~ 5i{M{Q)). The invertibility of M{Q) 
together with the fact that it has a block LU decomposition of type d shows that 5i{M{Q)) is invertible in L 
and we are done. □ 

We are now in position to argue as above. For all m, write ujMm = LmUm the block LU decomposition 
of Mm of type d. By the Chinese Remainder Theorem, there exists a unit lower triangular matrix L with 
coefficients in A such that L = Lm (mod Xm'"''^^) for all m. The matrix M — uj^^L then satisfies the 
two assumptions of Lemma l2.18l Hence, it satisfies also the conclusions of this Lemma and we have solved 
our problem. Algorithm |5] summarizes the different steps of the proposed solution. Of course, if uj is not 
invertible or one of the cjAfm's does not admit a block LU decomposition of the required type, Algorithm|5] 
fails. If it happens, we simply rerun the algorithm again and again until it works: it follows from Theorem 
I2.12l that we will get the desired answer quite fast. 

Let us analyze quickly how much precision is loss in average by this method. In order to fix ideas, let 
us assume that the entries of the matrix Af,„ are explicitely given as polynomials in Km[Xm] (eventually 
modulo ''^ ) and that all these polynomials are known with precision 0{tt^) for some integer N. For 
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Algorithm 5: A solution to Problem l2.17l 



1 D ^ Dmg{Pl' ' ■ ■ ■ . . , • • • Pn"'"); 

2 u <— a random matrix in Md{R); 

3 for mfrom 1 to n do 

4 compute d„ = {dm,l, ■ • ■ , dm,r^)\ 

s Lm <— Ld {ioMm) (computed by Algorithm |4|l; 

6 L a unit lower triangular matrix in Md{A) such that L = L,„ (mod Xm" '*^^) for all to; 

7 return {uj^^L, D); 



simplicity, we assume moreover that Mm(0) has coefficients in the ring of integers R„i of Km and that it 
is invertible in Md{Rmj^ Set as before jm{s) — rfm.i + . . . + c?m,s and, for all admissible pair (m, s), let 
denote by D,rL,s the determinant of the j„i(s)-th principal minor of a;A/„i. Define also 5{ljj) to the maximum 
of all fi?(-Dm,s(0)) when m and s run over all the possibilities. By the proof of Theorem 12. 121 we know 
that 5{ijj) is less than 

^ = logq(n^) + logg(''i + ---+rn) 

with probability at least i. In many concrete situations, it is not easy to compute exactly the r„j's but it 
will nevertheless in general quite simple to estimate them. Indeed, going back to the definition, it is clear 
that is less than both d and e.m,d and these latter quantites are natural parameters on which we will in 
general have a good control (cf |i2J, §3.2 for a concrete example). From now on, we assume that all matrices 
M„i computed on Une |5] satisfy this estimation. If this property does not hold, we simply agree to rerun 
Algorithm|5]until the desired property holds. 

The next step is to measure the size of the denominators appearing in the following nonconstant co- 
efficients. In order to do this, we introduce a new parameter w by requiring that all matrices Mm have 
coefficients in the ring Rm.iu defined as the image of -Rm[^%] in the quotient ring Km[Xm\/ X^'"^^^ . 
Clearly Dm.s belongs to Rm,w for all (to, s) and, by we have said before, it has a representant whose con- 
stant coefficient has a valuation less than v. Its inverse D^^^^ then belongs to tt^^ • Rm.v+w and is known 
up to an element of tt^^^" • Rm.v+w All entries of L,„ will consequently be known with this precision. 

It remains to analyze the line |6] of Algorithm |5] Note that the matrix L we want to compute can be 
expressed in terms of the L,„'s by the formula L = CiL[ + ■ • • + CnL'^ where: 

• L'm is a matrix with coefficients in K[X] whose reduction modulo Pm'™^^ corresponds to Lm via 
the natural isomorphism: 

K[X]/P:^'-+' ^ Km[Xm]lX'^-+\ Fix) ^ V E^^lpA . XI (13) 

i=0 

• Cm is a polynomial congruent to 1 modulo Pm'"^^^ and divisible by i^^jl '"^^ for all to' ^ m. 

In order to bound the loss of precision as we would like to do, we assume for simplicity that all Pm's 
are entirely known. We introduce again two new parameters. The first one is an integer vi for which we 
require that the image of R[X]/Pm''"^^ under the isomorphism ( fT3] ) contains n""^ ■ ■ ■ Rm[Xm]/ Xm'"^^^ 
for all to. The second parameter is the integer V2 defined as the opposite of the smallest valuation of 
a coefficient of the unique polynomial Cm of degree < J2m=ii^d.jn + l)degPm satisfying the above 
condition. Now, remember that we have proved that Lm are known up to an element of tt^"^" • Rm.v+w It 
is then a fortiori known up to an element of T^N-2v-e(v+w) ^vjjgj-g g = max(ei_d, • • ■ , e„,d). Inverting the 
isomorphism ( flj] ). we find that L'm is certainly known modulo 7r^~2ii-e„,d(f+«')-i'i Finally the formula 
L = CiL[ + ... + C„L; shows that L is known with precision (9(7r^-2t,-e(t,+u.)-ui-i,2) (j-g^-^^jj jjjat we 
have assumed that the P,„'s — and consequently the C„i's — are known with infinite precision). The total 
loss of precision of Algorithm|5]is then bounded by 2v + e{v + w) + vi + V2- 

Remark 2.20. The parameters vi and V2 are not easy to estimate in general. One can nevertheless keep in 
mind the following: vi measures the ramification of the roots of the Pm's and measures the distance 

^Otherwise, we would need to take in account the valuation of det Mm{0) as in Remarks l2.lQl and l2.13l 
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between these roots. For instance, to be more precise, one can easily prove that if all am lie in the ring 
of integers of an unramified extension K' of K then one can just take vi = 0. If in addition the a„i's are 
pairwise distinct in the residue field of K' (which is the same as to be distinct in the residue field of K since 
K' / K is unramified), one can also take V2 = 0. In that very particular case, the computation of line|6]does 
not generate any loss of precision. We refer to JS] for a quite different example where the constants vi and 
V2 do not vanish but stay nevertheless under control. 

Here is a final important remark. Algorithm |5] still works if, instead of computing (the i-part) of the 
block LU decomposition of LuMm, we compute a unit lower triangular (and not block unit lower triangular) 
Lm such that there exists a block upper triangular (with respect to d) matrix Um with the property that 
Mm = LmUm- Indeed, the knowledge of these L„i's is enough to compute L (which need to be only unit 
lower triangular) and then to conclude using Lemma 12.181 This remark is important because Algorithm!!] 
spends some time in line [TO] in clearing entries in order to make the computed matrix L block unit lower 
triangular instead of simply unit lower triangular. In other words, commenting the line [TO] in Algorithm!?] 
speeds up the execution of Algorithm!5]but do not have any influence on its correctness. 
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